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ABSTRACT 

Q ■ The gravitational-wave and accretion driven evolution of neutron stars in 

low mass X-ray binaries and similar systems is analyzed, while the amplitude 

c^ ■ of the radiating perturbation (here assumed to be an r-mode) remains small. If 

most of the star is superfluid, with (temperature independent) mutual friction 
^ ■ dominating the ordinary (temperature dependent) shear viscosity, the amplitude 

j^ ■ of the mode and the angular velocity of the star oscillate about their equilibrium 

values with a period of at least a few hundred years. The resulting oscillation 
of the neutron star temperature is also computed. For temperature dependent 
viscosity, the general conditions for the equilibrium to be stable are found. 
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Introduction 



We shall study the evolution of rapidly rotating accreting neutron stars under the influ- 
ence of their emission of gravitational radiation. We modify and extend the two-component 
model of the star (equilibrium plus perturbation) introduced by Owen et al. (1998) and 
also employed by Levin (1999), but restrict the analysis to small perturbations. These are 
assumed to be in the form of r-modes (Andersson 1998; Friedman & Morsink 1998; Lind- 
blom, Owen & Morsink 1998; Andersson, Kokkotas & Schutz 1999), which radiate mainly 
via Coriolis-driven velocity perturbations rather than the density perturbations of the less 
powerful f-modes. We must allow for large uncertainties in many of the relevant properties 
of neutron stars, such as the superfiuid transition temperature and the effects of magnetic 
fields and a core-crust boundary layer. 

After developing a general formalism, we shall study the evolution of the accreting 
neutron star under conditions in which its angular velocity remains approximately constant. 
One reason for this restriction is our interest in conditions in which the gravitational radiation 
is persistent over long time scales. This requires the existence of a stable or (possibly) 
overstable equilibrium. We shall see under what conditions the system can evolve toward 
this equilibrium, in which the rate of accretion of angular momentum from the surrounding 
disk is balanced by its rate of loss via gravitational radiation. If this equilibrium is achieved, 
the observed flux of gravitational radiation can be shown to be proportional to the observed 
flux of X-rays from the accretion (Wagoner 1984; Bildsten 1998). 

One of our longer term goals is the development of parameterized expressions describing 
possible time evolutions of the gravitational-wave frequency and amplitude, to facilitate 
detection by LIGO, VIRGO, and similar laser interferometer detectors. The brightest low 
mass X-ray binaries (LMXBs) are the prime targets. 



2. Dynamical and Thermal Evolution 

Consider a Newtonian neutron star in equilibrium (with equatorial radius R) which is 
perturbed by a nonaxisymmetric infinitesimal fluid displacement ^ = /(r, ^)e*'^™'^'*''^*-' ~ aR, 
with a ^ 1. Based on the work of Friedman & Schutz (1978a) and Levin & Ushomirsky 
(2001a), the total angular momentum J of the star can be decomposed into its equilibrium 
angular momentum J* and a perturbation proportional to the canonical angular momentum 
Jc- That is, 

J = J,iM, fi) + (1 - Kj)J, , J, = -K.a^J, , (1) 
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where M is the mass and Vt is the (uniform) angular velocity of the equilibrium star. All 
constants K(^ ) will be dimensionless, with Kj r^ K^. ^ 1. 

In classical mechanics, the action / = E /u of any normal mode of a set of oscillators 
(with frequency u) is an adiabatic invariant. For a fluid, the analogous quantity should 
be Ec/oo, where Ec is the canonical energy of the perturbation in the corotating frame and 
u = a + mVt is its frequency in that frame. However, we also have the general relation 
Ec = —{uj/m)Jc (Friedman & Schutz 1978a). Therefore, following Ho & Lai (2000), we 
assume that the canonical angular momentum is also an adiabatic invariant, and should 
therefore be unaffected by the slow rate of mass accretion. Thus it obeys the usual relation 
(Friedman & Schutz 1978b) 

dJjdt = 2J,[(F,(M, n) - F,{M, n, T„)] , (2) 

where Eg is the gravitational radiation growth rate and F„ is the viscous damping rate. The 
latter usually depends upon a spatially averaged temperature T^(t). 

On the other hand, conservation of total angular momentum requires that 

dJ/dt = 2J,Eg + Jait) , (3) 

where Ja = jaM is the rate of accretion of angular momentum. The mass is accreted with 
specific angular momentum ja at a rate M{t). 

Combining these equations then gives the dynamical evolution relations 

l^ = Eg-E^+[K,Eg + {l-K,)E^]Kca'-(^^^M{t), (4) 

{ja-j*y 



jj^ = -2[K,Eg + {l-K,)E,]Kca' + 
where h{M, il) = dJjdVt and j;(M, Vt) = dJjdM. 



'J *: 



M{t) ; (5) 



In obtaining our thermal evolution relation, the large uncertainties in some thermo- 
dynamic properties of the neutron star make it sufficient to consider slowly rotating stars. 
Thermal energy conservation for the entire star then gives 

-g^cJV = C{T)— = 2E,E,{T,) + Kn{M)c^ - L,{T,) , (6) 

where the rotating frame canonical energy E^ = KS^J^a^. Since the main contributor to 
the specific heat is the degenerate relativistic electrons, its value at constant volume (c^,) is 
essentially the same as that at constant pressure. The emissivities are produced by viscous 
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heating, pycnonuclear reactions and neutron emissions in the inner crust (proportional to a 
time-averaged mass accretion rate), and neutrino losses. The hydrogen/helium burning rate 
is assumed to be balanced by the surface emission of photons (Schatz et al. 1999), especially 
at the large accretion rate M = 10~^Mq yr~^ that we shall use. The mass accretion rate 
can be estimated from accretion energy conservation. The photon luminosity arising directly 
from the accretion is Lace ~ {GM/R)M{t), for a slowly rotating neutron star with a negligible 
magnetosphere. 

We are interested in the evolution of neutron stars after they have been spun up to the 
point where the gravitational radiation growth rate has become equal to the viscous damping 
rate: 

F,(fio, Mo) = F,(fio, Mo, To) = Fo , (7) 

so the evolution equation (2) vanishes. This equality defines our initial state. Before that 
time, we see from equation(2) that any intrinsic perturbation could not grow from its (in- 
finitesimal) value amin- The initial temperature To is determined by the vanishing of equation 
(6), with the nuclear heating in the inner crust balanced by the neutrino emission (Brown 
2000). Since we are only considering conditions in which a^ <^ 1, the properties Q{t) and 
M{t) evolve much slower than a{t) and T(t). For any other property Q* of the unperturbed 
star, let Qo = <5*(Mo,i7o). 

From now on we shall take the perturbation to be due to the dominant I = m = 2 
r-mode. In order to facilitate comparison with previous results, we shall adopt the neutron 
star model of Owen et al. (1998) for numerical work. Then the gravitational radiation growth 
rate of this mode is 

Fg = n^Tgr , V = 3-26 sec , n = VI{'kG{p)Y^''^ . (8) 

In the temperature range of interest (10^ < T < 10^° K), the viscous damping rate of this 
mode is approximated as 



'mf 



T,= — e-(^^/^^) +- ^^-^ , (9) 



Tsh V T, 



where Tc is the superfluid transition temperature. 

The first term is the contribution from the mutual friction between the neutron su- 
perfluid and the superconducting proton-relativistic electron fluid (Lindblom & Mendell 
2000). Its behavior (and that of other properties considered below) as the temperature 
passes through the superfluid transition temperature Tc is approximated by the exponential. 
Although Lindblom & Mendell (2000) found that Tmf ^ 10"^ sec, we shall keep it as a free pa- 
rameter because of the many uncertainties involved in the physics of this system, especially 
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when including magnetic effects such as the interaction of the vortex hues and flux tubes 
(Ruderman 2000). In fact, we fix it to satisfy our initial condition (7), when observationally 
relevant values of Qq, Mq, Tq are chosen. To approximately match an inferrred maximum spin 
rate of 330 Hz for the neutron stars in LMXBs, we shall choose CIq = VLQ{nG{p))~^/'^ = 0.25 
(Levin 1999), which then fixes Tmf ~ 13 sec if T <^ Tc. The maximum rotation rate of a 
neutron star corresponds to fi = 2/3. 

The second term is the contribution of the ordinary shear viscosity. Lindblom & Mendell 
(2000) obtained Tsh = 1.0 x 10^ sec for their superfluid neutron star model, dominated by 
electron-electron scattering in the superfluid regions and neutron-neutron scattering in the 
normal regions. Above the superfluid transition temperature, Lindblom, Owen & Morsink 
(1998) obtained Tgh = 2.5 x 10® sec. There are also large uncertainties in this contribution, due 
to the shear in a boundary layer between the core and crust (Andersson, Jones, Kokkotas 
& Stergioulas 1999; Bildsten & Ushomirsky 2000; Lindblom, Owen & Ushomirsky 2000; 
Rieutord 2000; Wu, Matzner & Arras 2001) and the uncertain response of the crust to the 
mode (Levin & Ushomirsky 2001b). 

In keeping with the fact that we are only working to lowest order in Cl (as well as the 
relativity parameter GM/Rc^), we take ja — Jo — ja and Jq = /q^o- We also note that 
Kc = 3Ke = 0.094. (The value of Kj is unimportant.) 

Now that we have specified all properties in the equations (4) and (5) of evolution of 
a{t) and fi(t), we can consider the thermal evolution. In what follows we shall assume that 
thermal conductivity timescales are short enough to give relations T^{T) and T,y{T) between 
these three spatially averaged temperatures that appear in equation (6). The normal and 
superfluid contribution to the specific heat and the neutrino luminosity that appear in this 
equation are approximated by 

^ V-' / ~ V-'norm^ "r^superJ-' K^norm ~ '^^^ super ) i \^^) 

Li^\T) = Ljjj^(j/r e " + L^^^^T . (11) 

The constants L'^jiqa ^"^^ ^'brem ^^^ obtained by fitting the results of Brown (2000) for normal 
and superfluid neutron stars. We also take the nuclear heating constant Kn = 1 x 10"^ 
(Brown 2000). 

In Figures 1 and 2 we show the results of integrating our three coupled evolution equa- 
tions (4), (5), and (6), if Tc ^ Tq. The main feature is the spin-down (due to loss of angular 
momentum in gravitational waves) and viscous heating during the time when the amplitude 
a of the mode exceeds its equilibrium value (denoted by the dashed line). 

For such large values of T^., the evolution qualitatively resembles that which occurs 
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Fig. 1. — The evolution of the mode amphtude a, for the case T^^ Tq. The dashed hne is 
its equihbrium value. Time is in units of 10^ years. Subsequent cycles are similar. 
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Fig. 2. — Same as Figure 1, for the angular velocity Q{nG{p)) ^^^ and the temperature 
T/10^ K. Here the dashed line is the initial value. 



when dF^/dT = (requiring neglect of the shear viscosity). In this case equations (4) and 
(5) decouple from equation (6). In addition, a; = Ina completely decouples, obeying the 
equation 



d'x/dr - 2K^Foe'''dx/dt + dV/dx = , V{x) ^ Fo{K,Foe''' - F^x 



(12) 



where Fa is the time-averaged rate of accretion of angular momentum. The Eddington limit 
to the accretion rate gives r^ = XjFa ^ 5 x 10® years. The sign of the damping term 
is opposite to that obtained by Levin (1999), leading to overstable oscillations about the 
equilibrium amphtude 

ae = (Fa/2ircFo)'/'<2xlO-^ 

with a period 

P ^ [8\n{ama./armn)/{FoFa)]'^^ > 300 years . 
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The fraction F> of the time that a > ae is 

F> ^ \n[8ln{a„,ax/amin)]/4:\n{amax/amin) ~ 0.1 . 

It is also found that amax increases and amin decreases on a timescale ~ P/amax- However, 
O-min is presumably limited by the intrinsic fluctuations in the neutron star, which is also why 
we have not allowed a to drop below its initial value in Figure 1. Because of the inclusion 
of the shear viscosity, the value of the period P in Figure 1 is about ten times longer than 
the above estimate while the value of the fraction F> in Figure 1 is about ten times less. 

On the other hand, for values of Tc < Tq, the temperature dependence of the viscous 
damping rate F„ is strong enough to produce a thermal runaway to large values of T and 
a (outside the range of validity of our equations), as found by Levin (1999). This result is 
generalized in the next section. 



3. Behavior Near Equilibrium 

In contrast to the initial state, the equilibrium state X* of our dynamical variables 

x\t) = {a,n,T} = xi[i + c{t)] , ICi « 1 , 

is defined by the vanishing of the evolution equation (3), in addition to equations (2) and 
(6). Employing a constant (averaged) accretion rate, the evolution equations give 



dC/dt = A'^C\ Cocexp(At), \\A'^ - X6'^\\ = . 
Assume now that |(9F„/(9T| ~ F^/Te, etc. 

The coefficients of the eigenvalue equation A^ + a2A^ + aiA + ao = are 



a2 



Oi 



ao 



1 



(§1 -(-. 



or 



KralFo , 






4K,n,alFo 

a. 



d{F, - F,) 

on 



[df 



16K,{E,),alF^ fdF, 



a 



dT 



~ KratFi^ 



(13) 



(14) 
(15) 



(16) 



The ratio of rotational to thermal energy is Kr = 2Ke^eJo/CeTe ~ 10^. We have used the 
fact that the coohng rate Fc = Ly(Tf.)/CeTe ~ KrolFQ. 



Now we also employ the inequalities Kr ^ 1 and (mode energy) /(thermal energy) 
~ Krai < 10""^ ^ |ai| ^ a^ to obtain the eigenvalues 

Ai,2 = -02/2 ± V-«i , A3 = -ao/oi . 

We have used the fact that IA3I ~ a^Fo <^ |Ai| ~ IA2I. 

Let us examine the two relevant possibilities. For cases such as dominance by shear 
viscosity in equation (9), with Ec> 0, 

a, oc {E,),{dFjdT), < ^ Ai,2 = ±V^ ~ K'^a^F^ • 

Thus this equilibrium is unstable, with a growth rate Ai that is of the same magnitude as 
found by Levin (1999). 

The other possibility is 

fli > =^ Ai^2 — —02/2 ± i^/oi . 

Thus stability also requires that ao > and 02 > 0. From their relations above, we see that 
this means that the variation in the cooling rate with temperature must be greater than 
twice the variation in the viscous heating rate with temperature (which is usually satisfied). 

This work was supported in part by NSF grant PHY-0070935 to R.V.W. and NASA 
grant NAS 8-39225 to Gravity Probe B. R.V.W thanks the Aspen Center for Physics for 
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2000 program on Spin and Magnetism in Young Neutron Stars at the Institute for Theoretcial 
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